!!!! Business_Cycles_Insurance Toolbox - (12/21/2022) !!!!
!! Subroutines for "The Impact of Private and Public Insurance on Consumption Over the Business Cycle"

!!!! Table of Contents !!!!
!! 1.) Forecasting Rules
!! 2.) Aggregate Fluctuations
!! 3.) UI Polices
!! 4.) Employment Transition Matrix
!! 5.) Individual State Transition Matrix
!! 6.) Aggregate State Transition Matrix
!! 7.) Utility Function (CRRA)
!! 8.) Borrowing Limit
!! 9.) Grid Search Subroutine

module BCI_TB

use BCI_Par
use Fortran_TB

implicit none

    contains

!!Introduction to Code
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 3.) UI Policies

subroutine UI_policies(switch_scenario, upsilon_d, upsilon_r, upsilon_bar)

!!!! Variables and Arrays !!!!
integer, intent(in) :: switch_scenario
real(rp), dimension(Nz), intent(out) :: upsilon_d, upsilon_r, upsilon_bar

!!!! UI Policies !!!!
upsilon_d(1) = upsilon_dg
upsilon_r(1) = upsilon_rg
upsilon_bar(1) = upsilon_bar_g
if ( switch_scenario <= 4 .or. switch_scenario == 8 .or. switch_scenario == 12 .or. switch_scenario == 14) then                 !! UI Extension
    upsilon_d(2) = upsilon_db
    upsilon_r(2) = upsilon_rg
    upsilon_bar(2) = upsilon_bar_g
elseif ( switch_scenario == 5 .or. switch_scenario == 9 .or. switch_scenario == 13 .or. switch_scenario == 15 ) then            !! Acyclical UI
    upsilon_d(2) = upsilon_dg
    upsilon_r(2) = upsilon_rg
    upsilon_bar(2) = upsilon_bar_g
elseif (switch_scenario == 6 .or. switch_scenario == 10 ) then                                                                  !! UI RR Increase
    upsilon_d(2) = upsilon_dg
    upsilon_r(2) = 0.62_rp
    upsilon_bar(2) = upsilon_bar_g
elseif ( switch_scenario == 7 .or. switch_scenario == 11 .or. switch_scenario == 16 .or. switch_scenario == 17 ) then           !! No UI
    upsilon_d(2) = upsilon_db
    upsilon_r(:) = zero
    upsilon_bar(2)= upsilon_bar_b
elseif (switch_scenario == 12 ) then                                                                                            !! UI Level Increase
    upsilon_d(2) = upsilon_dg
    upsilon_r(2) = upsilon_rb
    upsilon_bar(2) = upsilon_bar_b
end if

end subroutine UI_policies

!! 3.) UI Policies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 4.) Employment Transition Matrix

subroutine emp_transitions(n, xi, lambda, upsilon_d, pi)

!!!! Variables and Arrays !!!!
integer, intent(in) :: n
real(rp), intent(in) :: xi, lambda, upsilon_d
real(rp), dimension(n,n), intent(out) :: pi

!!!! Transition Matrix !!!!
pi = zero
pi(1,1) = one - xi
pi(1,2) = xi
pi(1,3) = zero
pi(2,1) = lambda
pi(2,2) = (one - upsilon_d)*(one - lambda)
pi(2,3) = upsilon_d*(one - lambda)
pi(3,1) = lambda
pi(3,2) = zero
pi(3,3) = one - lambda

end subroutine emp_transitions

!! 4.) Employment Transition Matrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 5.) Individual Transition Matrix

subroutine ind_transitions(ni, ne, nn, nb, pi_e, pi_n, e_index, n_index, b_index, i_index, pi_i )

!!!! Variables and Arrays !!!!
integer, intent(in) :: ni, ne, nn, nb
real(rp), dimension(ne,ne), intent(in) :: pi_e
real(rp), dimension(nn,nn), intent(in) :: pi_n
integer :: i, e, n, b, i2, e2, n2, b2, count
integer, dimension(ni), intent(out) :: e_index, n_index, b_index
integer, dimension(ne,nn,nb), intent(out) :: i_index
real(rp), dimension(ni,ni,2), intent(out) :: pi_i

!!!! Assign Indices !!!!
count = 0
do e = 1,ne
	do n = 1,nn
		do b = 1,nb
			count = count + 1
			e_index(count) = e
			n_index(count) = n
			b_index(count) = b
            i_index(e,n,b) = count
		end do
	end do
end do

!!!! Overall Transition Matrix !!!!
pi_i = zero
do i = 1,ni
	n = n_index(i)
	e = e_index(i)
	b = b_index(i)

	do i2 = 1,ni
		n2 = n_index(i2)
		e2 = e_index(i2)
		b2 = b_index(i2)
        
		if (b == b2)  then
			if ( n == 1 ) then
				if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 2 .and. e == e2 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                elseif (n2 == 3 .and. e == e2 ) then
                    pi_i(i,i2,2) = one
                endif
            elseif (n == 2 ) then
                if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 2 .and. e == e2) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                elseif ( n2 == 3 .and. e == e2 )  then
                    pi_i(i,i2,1) = pi_n(n,n2)
                    pi_i(i,i2,2) = one
                end if
            else
                if ( n2 == 1 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)*pi_e(e,e2)
                elseif ( n2 == 3 .and. e == e2 ) then
                    pi_i(i,i2,1) = pi_n(n,n2)
                    pi_i(i,i2,2) = one
                end if
            end if
		end if
        
	end do
end do

end subroutine ind_transitions

!! 5.) Individual Transition Matrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 6.) Aggregate Transition Matrix

subroutine agg_transitions( nz, rho_z1, rho_z2, pi_z, cum_z )

!!!! Variables and Arrays !!!!
integer, intent(in) :: nz
real(rp), intent(in) :: rho_z1, rho_z2
real(rp), dimension(nz,nz), intent(out) :: pi_z, cum_z

!!!! Aggregate Transition Matrix !!!!
pi_z = zero
cum_z = zero
pi_z(1,1) = rho_z1
pi_z(1,2) = one - rho_z1
pi_z(2,1) = one - rho_z2
pi_z(2,2) = rho_z2
cum_z(1,1) = pi_z(1,1)
cum_z(1,2) = pi_z(1,1) + pi_z(1,2)
cum_z(2,1) = pi_z(2,1)
cum_z(2,2) = pi_z(2,1) + pi_z(2,2)

end subroutine agg_transitions

!! 6.) Aggregate Transitions
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 7.) Utility Function

subroutine utility(beta, chi, consumption, EV, u_val)

!!!! Variables and Arrays !!!!
real(rp), intent(in) :: beta, chi, consumption, EV
real(rp), intent(out) :: u_val

!!!! Function !!!!
u_val = (one/(one-sigma))*(( consumption )**(one-sigma)) - chi + beta*EV

end subroutine utility

!! 7.) Utility Function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! 8.) Borrowing Limit

subroutine borrowing_limit(nd, xgrid, discount, iB, minB)

!!!! Variables and Arrays !!!!
integer, intent(in) :: nd
real(rp), dimension(nd), intent(in) :: xgrid, discount
integer, intent(out) :: iB
real(rp), intent(out) :: minB
integer :: x

!!!! Calculate Endogenous Borrowing Limit !!!!
minB = 0.0_rp
iB = nd
do x = 1,nd
	if ( xgrid(x)*discount(x) < minB ) then
		minB = xgrid(x)*discount(x)
		iB = x
	end if
end do

end subroutine borrowing_limit

!! 8.) Borrowing Limit
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 9.) Grid Search

subroutine grid_search(ibar, beta, chi, a_val, income, discount, agrid, EV, maxV, maxG, maxC, maxQ)

!!!! Variables and Arrays !!!!
integer, intent(in) :: ibar
real(rp), intent(in) :: beta, chi, a_val, income
real(rp), dimension(nd), intent(in) :: discount
real(rp), dimension(na), intent(in) :: agrid, EV
real(rp), intent(out) :: maxV, maxG, maxC, maxQ
integer :: a2
real(rp) :: consumption, Value

!!!! Optimal Value for Borrowing !!!!
maxV = disutility; maxG = zero; maxC = zero; maxQ = zero
do a2 = ibar,nd
    consumption = a_val + income + y_bar - discount(a2)*agrid(a2)
    if ( consumption > zero ) then
        Value = (one/(one-sigma))*(( consumption )**(one-sigma)) - chi + beta*EV(a2)
    else
        Value = disutility
    end if

    if ( Value > maxV ) then
    	maxV = Value
    	maxG = agrid(a2)
    	maxC = consumption
        maxQ = discount(a2)
    end if
end do

end subroutine grid_search

!! 9.) Grid Search
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 10.) Income Distribution

subroutine income_distribution(nx, e_index, n_index, upsilon_r, gamma_j, egrid, i_income, j_income, income_grid)

!!!! Variables and Arrays !!!!
integer :: i, i2, j, e, n, count, ix
real(rp) :: income, omega_x
integer, intent(in) :: nx       !!dimension of income grid (ni*nj)
integer, dimension(ni), intent(in) :: e_index, n_index
real(rp), intent(in) :: upsilon_r
real(rp), dimension(nj), intent(in) :: gamma_j
real(rp), dimension(ne), intent(in) :: egrid
integer, dimension(nx), intent(out) :: i_income, j_income
real(rp), dimension(nx) :: income_grid

!!!! Initialize Income Grid !!!!
i_income(:) = 1
j_income(:) = 1
income_grid(:) = zero

!!!! Income Order - Positive Income !!!!
count = 1
do j = 1,nj
    do i = 1,ni
        e = e_index(i); n = n_index(i)
        
        if ( n == 1 ) then
            income = gamma_j(j)*egrid(e)
        elseif ( n == 2 ) then
            income = upsilon_r*gamma_j(j)*egrid(e)
        else
            income = zero
        end if
        
        if ( count < 3 ) then
            i_income(count) = i
            j_income(count) = j
            income_grid(count) = income
            count = count + 1
        else
            
            if ( income >= income_grid(count-1) ) then
                i_income(count) = i
                j_income(count) = j
                income_grid(count) = income
                count = count + 1
            elseif ( income <= income_grid(1) ) then
                do i2 = count,2,-1
                    i_income(i2) = i_income(i2-1)
                    j_income(i2) = j_income(i2-1)
                    income_grid(i2) = income_grid(i2-1)
                end do
                i_income(1) = i
                j_income(1) = j
                income_grid(1) = income
                count = count + 1
            else
                call interpolation(income_grid, income, 1, count-1, nx, ix, omega_x)
                do i2 = count,ix+2,-1
                    income_grid(i2) = income_grid(i2-1)
                    i_income(i2) = i_income(i2-1)
                    j_income(i2) = j_income(i2-1)
                end do
                income_grid(ix+1) = income
                i_income(ix+1) = i
                j_income(ix+1) = j
                count = count + 1
            end if
            
        end if
        
    end do
end do


end subroutine income_distribution

!! 10.) Income Distribution
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!! 11.) Future Index - For Simulation

subroutine future_index(cum_pi, rand_x, nx, ix)

!!!! Variables and Arrays !!!!
integer :: x
real(rp) :: omega
integer, intent(in) :: nx
real(rp), intent(in) :: rand_x
real(rp), dimension(nx), intent(in) :: cum_pi
integer, intent(out) :: ix

!!!! Calculate Future Index !!!!
if ( rand_x <= cum_pi(1) ) then
    ix = 1
else
    call interpolation(cum_pi, rand_x, 1, nx, nx, x, omega)
    ix = x + 1
end if

end subroutine


end module BCI_TB